{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mendelian simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_sims = 100000\n",
    "num_ofs = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAIMCAYAAAAq64s8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X2MXXd95/HPd+OC+sQSGpNNk7BOkakEqE3BCqwQVXZTgkkqEqrSOrsqXsrKhCWrot2VMK20QXSR0gdaLSs2VSgWYQUJaWmKpYQGN1sVrURoHJrmAUjjpC4xsRKXsDyIiirhu3/MMVwmM+PBM/bMz/N6SVf33t85585vfHTH9z333DPV3QEAAIAR/bO1ngAAAAAcL1ELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwNq31BI7XGWec0Vu2bFnraQAAAHAC3HXXXf/Q3ZuPtd6wUbtly5bs379/racBAADACVBVf7+c9Rx+DAAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADCsY0ZtVe2pqser6r6ZsY9W1d3T5WBV3T2Nb6mqf5xZ9gcz27y0qu6tqgNV9d6qqmn8OVW1r6oenK5PPxHfKAAAAKee5bxT+8Ek22cHuvuXu/v87j4/yceS/MnM4oeOLuvuK2fGr02yK8nW6XL0MXcnub27tya5fboPAAAAx3TMqO3uTyV5YqFl07utv5TkhqUeo6rOSvKs7v50d3eSDyW5fFp8WZLrp9vXz4wDAADAklb6mdpXJnmsux+cGTuvqv66qv6yql45jZ2d5NDMOoemsSQ5s7sPJ8l0/dwVzgkAAIANYtMKt78i3/su7eEkz+vuL1fVS5P8aVW9KEktsG1/v1+sqnZl7hDmPO95zzuO6QIAAHAqOe6orapNSX4hyUuPjnX3t5J8a7p9V1U9lOQFmXtn9pyZzc9J8uh0+7GqOqu7D0+HKT++2Nfs7uuSXJck27Zt+76jGAD4Xlt237LWU1iWg9dcutZTAGCdWsnhxz+X5Avd/Z3Diqtqc1WdNt3+icydEOrh6bDir1fVy6fP4b4hycenzfYm2Tnd3jkzDgAAAEtazp/0uSHJp5P8ZFUdqqo3TYt25OkniPrZJPdU1d8k+eMkV3b30ZNMvSXJHyY5kOShJJ+Yxq9J8qqqejDJq6b7AAAAcEzHPPy4u69YZPzfLzD2scz9iZ+F1t+f5MULjH85yUXHmgcAAADMt9KzHwMAAMCaEbUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxr01pPAABONVt237LWUzjljPRvevCaS9d6CgAbindqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYR0zaqtqT1U9XlX3zYy9s6q+VFV3T5dLZpa9o6oOVNUDVfXqmfHt09iBqto9M35eVX2mqh6sqo9W1TNW8xsEAADg1LWcd2o/mGT7AuO/393nT5dbk6SqXphkR5IXTdv8r6o6rapOS/K+JK9J8sIkV0zrJslvTY+1NclXkrxpJd8QAAAAG8cxo7a7P5XkiWU+3mVJbuzub3X33yU5kOSC6XKgux/u7n9KcmOSy6qqkvybJH88bX99ksu/z+8BAACADWoln6m9qqrumQ5PPn0aOzvJIzPrHJrGFhv/sST/r7ufnDcOAAAAx3S8UXttkucnOT/J4STvmcZrgXX7OMYXVFW7qmp/Ve0/cuTI9zdjAAAATjnHFbXd/Vh3P9Xd307y/swdXpzMvdN67syq5yR5dInxf0jy7KraNG98sa97XXdv6+5tmzdvPp6pAwAAcAo5rqitqrNm7r4uydEzI+9NsqOqnllV5yXZmuSvktyZZOt0puNnZO5kUnu7u5P8RZJfnLbfmeTjxzMnAAAANp5Nx1qhqm5IcmGSM6rqUJKrk1xYVedn7lDhg0nenCTdfX9V3ZTkc0meTPLW7n5qepyrktyW5LQke7r7/ulLvD3JjVX135P8dZIPrNp3BwAAwCntmFHb3VcsMLxoeHb3u5O8e4HxW5PcusD4w/nu4csAAACwbCs5+zEAAACsKVELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwjhm1VbWnqh6vqvtmxn6nqr5QVfdU1c1V9expfEtV/WNV3T1d/mBmm5dW1b1VdaCq3ltVNY0/p6r2VdWD0/XpJ+IbBQAA4NSznHdqP5hk+7yxfUle3N0/leRvk7xjZtlD3X3+dLlyZvzaJLuSbJ0uRx9zd5Lbu3trktun+wAAAHBMx4za7v5UkifmjX2yu5+c7t6R5JylHqOqzkryrO7+dHd3kg8luXxafFmS66fb18+MAwAAwJJW4zO1v5rkEzP3z6uqv66qv6yqV05jZyc5NLPOoWksSc7s7sNJMl0/d7EvVFW7qmp/Ve0/cuTIKkwdAACAka0oaqvqN5I8meTD09DhJM/r7p9J8p+TfKSqnpWkFti8v9+v193Xdfe27t62efPm4502AAAAp4hNx7thVe1M8vNJLpoOKU53fyvJt6bbd1XVQ0lekLl3ZmcPUT4nyaPT7ceq6qzuPjwdpvz48c4JAACAjeW43qmtqu1J3p7ktd39zZnxzVV12nT7JzJ3QqiHp8OKv15VL5/OevyGJB+fNtubZOd0e+fMOAAAACzpmO/UVtUNSS5MckZVHUpydebOdvzMJPumv8xzx3Sm459N8q6qejLJU0mu7O6jJ5l6S+bOpPyDmfsM7tHP4V6T5KaqelOSLyZ5/ap8ZwAAAJzyjhm13X3FAsMfWGTdjyX52CLL9id58QLjX05y0bHmAQAAAPOtxtmPAQAAYE2IWgAAAIYlagEAABiWqAUAAGBYohYAAIBhHfPsxwAALN+W3bes9RSW7eA1l671FABWzDu1AAAADEvUAgAAMCxRCwAAwLBELQAAAMNyoigAhjDSyXcAgJPHO7UAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADCsZUVtVe2pqser6r6ZsedU1b6qenC6Pn0ar6p6b1UdqKp7quolM9vsnNZ/sKp2zoy/tKrunbZ5b1XVan6TAAAAnJqW+07tB5Nsnze2O8nt3b01ye3T/SR5TZKt02VXkmuTuQhOcnWSlyW5IMnVR0N4WmfXzHbzvxYAAAA8zbKitrs/leSJecOXJbl+un19kstnxj/Uc+5I8uyqOivJq5Ps6+4nuvsrSfYl2T4te1Z3f7q7O8mHZh4LAAAAFrWSz9Se2d2Hk2S6fu40fnaSR2bWOzSNLTV+aIFxAAAAWNKJOFHUQp+H7eMYf/oDV+2qqv1Vtf/IkSMrmCIAAACngpVE7WPTocOZrh+fxg8lOXdmvXOSPHqM8XMWGH+a7r6uu7d197bNmzevYOoAAACcClYStXuTHD2D8c4kH58Zf8N0FuSXJ/nqdHjybUkurqrTpxNEXZzktmnZ16vq5dNZj98w81gAAACwqE3LWamqbkhyYZIzqupQ5s5ifE2Sm6rqTUm+mOT10+q3JrkkyYEk30zyxiTp7ieq6jeT3Dmt967uPnryqbdk7gzLP5jkE9MFAAAAlrSsqO3uKxZZdNEC63aSty7yOHuS7FlgfH+SFy9nLgAAAHDUiThRFAAAAJwUohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGEdd9RW1U9W1d0zl69V1duq6p1V9aWZ8UtmtnlHVR2oqgeq6tUz49unsQNVtXul3xQAAAAbw6bj3bC7H0hyfpJU1WlJvpTk5iRvTPL73f27s+tX1QuT7EjyoiQ/nuTPq+oF0+L3JXlVkkNJ7qyqvd39ueOdGwAAABvDcUftPBcleai7/76qFlvnsiQ3dve3kvxdVR1IcsG07EB3P5wkVXXjtK6oBQAAYEmr9ZnaHUlumLl/VVXdU1V7qur0aezsJI/MrHNoGltsHAAAAJa04qitqmckeW2SP5qGrk3y/Mwdmnw4yXuOrrrA5r3E+EJfa1dV7a+q/UeOHFnRvAEAABjfarxT+5okn+3ux5Kkux/r7qe6+9tJ3p/vHmJ8KMm5M9udk+TRJcafpruv6+5t3b1t8+bNqzB1AAAARrYan6m9IjOHHlfVWd19eLr7uiT3Tbf3JvlIVf1e5k4UtTXJX2XundqtVXVe5k42tSPJv12FeQFwDFt237LWUwAAWJEVRW1V/VDmzlr85pnh366q8zN3CPHBo8u6+/6quilzJ4B6Mslbu/up6XGuSnJbktOS7Onu+1cyLwAAjm2UX2wdvObStZ4CsI6tKGq7+5tJfmze2K8ssf67k7x7gfFbk9y6krkAAACw8azW2Y8BAADgpBO1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMa8VRW1UHq+reqrq7qvZPY8+pqn1V9eB0ffo0XlX13qo6UFX3VNVLZh5n57T+g1W1c6XzAgAA4NS3Wu/U/uvuPr+7t033dye5vbu3Jrl9up8kr0mydbrsSnJtMhfBSa5O8rIkFyS5+mgIAwAAwGJO1OHHlyW5frp9fZLLZ8Y/1HPuSPLsqjoryauT7OvuJ7r7K0n2Jdl+guYGAADAKWI1oraTfLKq7qqqXdPYmd19OEmm6+dO42cneWRm20PT2GLj36OqdlXV/qraf+TIkVWYOgAAACPbtAqP8YrufrSqnptkX1V9YYl1a4GxXmL8ewe6r0tyXZJs27btacsBAADYWFb8Tm13PzpdP57k5sx9Jvax6bDiTNePT6sfSnLuzObnJHl0iXEAAABY1Iqitqp+uKp+9OjtJBcnuS/J3iRHz2C8M8nHp9t7k7xhOgvyy5N8dTo8+bYkF1fV6dMJoi6exgAAAGBRKz38+MwkN1fV0cf6SHf/WVXdmeSmqnpTki8mef20/q1JLklyIMk3k7wxSbr7iar6zSR3Tuu9q7ufWOHcAAAAOMWtKGq7++EkP73A+JeTXLTAeCd56yKPtSfJnpXMBwAAgI3lRP1JHwAAADjhRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMa9NaTwAAAJayZfctaz2FZTt4zaVrPQXYcEQtwAkw0gswAICROfwYAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFjHHbVVdW5V/UVVfb6q7q+qX5vG31lVX6qqu6fLJTPbvKOqDlTVA1X16pnx7dPYgaravbJvCQAAgI1i0wq2fTLJf+nuz1bVjya5q6r2Tct+v7t/d3blqnphkh1JXpTkx5P8eVW9YFr8viSvSnIoyZ1Vtbe7P7eCuQEAALABHHfUdvfhJIen21+vqs8nOXuJTS5LcmN3fyvJ31XVgSQXTMsOdPfDSVJVN07riloAAACWtCqfqa2qLUl+JslnpqGrquqeqtpTVadPY2cneWRms0PT2GLjAAAAsKQVR21V/UiSjyV5W3d/Lcm1SZ6f5PzMvZP7nqOrLrB5LzG+0NfaVVX7q2r/kSNHVjp1AAAABreiqK2qH8hc0H64u/8kSbr7se5+qru/neT9+e4hxoeSnDuz+TlJHl1i/Gm6+7ru3tbd2zZv3rySqQMAAHAKWMnZjyvJB5J8vrt/b2b8rJnVXpfkvun23iQ7quqZVXVekq1J/irJnUm2VtV5VfWMzJ1Mau/xzgsAAICNYyVnP35Fkl9Jcm9V3T2N/XqSK6rq/MwdQnwwyZuTpLvvr6qbMncCqCeTvLW7n0qSqroqyW1JTkuyp7vvX8G8AAAA2CBWcvbj/5uFPw976xLbvDvJuxcYv3Wp7QAAAGAhq3L2YwAAAFgLohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBYohYAAIBhiVoAAACGJWoBAAAYlqgFAABgWKIWAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGFtWusJACzXlt23rPUUAABYZ0QtAACskpF+AXvwmkvXegqwKhx+DAAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLBELQAAAMMStQAAAAxL1AIAADAsUQsAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADEvUAgAAMCxRCwAAwLA2rfUEAACAk2/L7lvWegrLdvCaS9d6CqxjohY2uJH+QwMAgPkcfgwAAMCwRC0AAADDErUAAAAMS9QCAAAwLFELAADAsEQtAAAAwxK1AAAADGvdRG1Vba+qB6rqQFXtXuv5AAAAsP5tWusJJElVnZbkfUleleRQkjuram93f25tZwbHb8vuW9Z6CgAAp4RRXlcdvObStZ7ChrRe3qm9IMmB7n64u/8pyY1JLlvjOQEAALDOrYt3apOcneSRmfuHkrxsjebCOjbKb+kAANh4Rnqteiq9q7xeorYWGOunrVS1K8mu6e43quqBEzqrlTkjyT+s9SRYFvtqHPbVOOyrcdhXY7CfxmFfjWND76v6rbWewbL8y+WstF6i9lCSc2fun5Pk0fkrdfd1Sa47WZNaiara393b1noeHJt9NQ77ahz21TjsqzHYT+Owr8ZhX5061stnau9MsrWqzquqZyTZkWTvGs8JAACAdW5dvFPb3U9W1VVJbktyWpI93X3/Gk8LAACAdW5dRG2SdPetSW5d63msoiEOkyaJfTUS+2oc9tU47Ksx2E/jsK/GYV+dIqr7aedjAgAAgCGsl8/UAgAAwPdN1K5QVW2vqgeq6kBV7V5g+TOr6qPT8s9U1ZaTP0uq6tyq+ouq+nxV3V9Vv7bAOhdW1Ver6u7p8t/WYq4kVXWwqu6d9sP+BZZXVb13el7dU1UvWYt5bnRV9ZMzz5e7q+prVfW2eet4Xq2RqtpTVY9X1X0zY8+pqn1V9eB0ffoi2+6c1nmwqnaevFlvPIvsp9+pqi9MP99urqpnL7Ltkj8rWV2L7Kt3VtWXZn7GXbLItku+XmR1LbKvPjqznw5W1d2LbOt5NSCHH69AVZ2W5G+TvCpzf5boziRXdPfnZtb5j0l+qruvrKodSV7X3b+8JhPewKrqrCRndfdnq+pHk9yV5PJ5++rCJP+1u39+jabJpKoOJtnW3Qv+7bjpRcN/SnJJkpcl+R/d/bKTN0Pmm34efinJy7r772fGL4zn1Zqoqp9N8o0kH+ruF09jv53kie6+ZnphfXp3v33eds9Jsj/Jtsz9zfi7kry0u79yUr+BDWKR/XRxkv8znUjzt5Jk/n6a1juYJX5WsroW2VfvTPKN7v7dJbY75utFVtdC+2re8vck+Wp3v2uBZQfjeTUc79SuzAVJDnT3w939T0luTHLZvHUuS3L9dPuPk1xUVXUS50iS7j7c3Z+dbn89yeeTnL22s2IFLsvcf1Td3Xckefb0iwvWzkVJHpoNWtZWd38qyRPzhmf/T7o+yeULbPrqJPu6+4kpZPcl2X7CJrrBLbSfuvuT3f3kdPeOJOec9InxNIs8p5ZjOa8XWUVL7avpdfgvJbnhpE6KE0rUrszZSR6ZuX8oTw+l76wz/Qf11SQ/dlJmx4KmQ8B/JslnFlj8r6rqb6rqE1X1opM6MWZ1kk9W1V1VtWuB5ct57nFy7cjiLxA8r9aPM7v7cDL3y74kz11gHc+v9eVXk3xikWXH+lnJyXHVdKj4nkUO6fecWl9emeSx7n5wkeWeVwMStSuz0Duu84/nXs46nCRV9SNJPpbkbd39tXmLP5vkX3b3Tyf5n0n+9GTPj+94RXe/JMlrkrx1OoxolufVOlJVz0jy2iR/tMBiz6vxeH6tE1X1G0meTPLhRVY51s9KTrxrkzw/yflJDid5zwLreE6tL1dk6XdpPa8GJGpX5lCSc2fun5Pk0cXWqapNSf55ju/QFVaoqn4gc0H74e7+k/nLu/tr3f2N6fatSX6gqs44ydMkSXc/Ol0/nuTmzB26NWs5zz1Ontck+Wx3PzZ/gefVuvPY0UP1p+vHF1jH82sdmE7Q9fNJ/l0vcgKUZfys5ATr7se6+6nu/naS92fhfeA5tU5Mr8V/IclHF1vH82pMonZl7kyytarOm96p2JFk77x19iY5eubIX8zciR/8du4kmz4/8YEkn+/u31tknX9x9PPOVXVB5p4fXz55syRJquqHp5N5pap+OMnFSe6bt9reJG+oOS/P3MkeDp/kqfJdi/7W2/Nq3Zn9P2lnko8vsM5tSS6uqtOnQykvnsY4Sapqe5K3J3ltd39zkXWW87OSE2ze+Rxel4X3wXJeL3Jy/FySL3T3oYUWel6Na9NaT2Bk01kJr8rcf/anJdnT3fdX1buS7O/uvZkLqf9dVQcy9w7tjrWb8Yb2iiS/kuTemVO4/3qS5yVJd/9B5n7p8JaqejLJPybZ4RcQa+LMJDdPHbQpyUe6+8+q6srkO/vq1syd+fhAkm8meeMazXXDq6ofytwZPd88Mza7rzyv1khV3ZDkwiRnVNWhJFcnuSbJTVX1piRfTPL6ad1tSa7s7v/Q3U9U1W9m7oV4krxY7K6kAAAAjklEQVSrux1hdIIssp/ekeSZSfZNPwvvmP6Kwo8n+cPuviSL/Kxcg29hw1hkX11YVedn7nDig5l+Fs7uq8VeL67Bt7BhLLSvuvsDWeD8D55XpwZ/0gcAAIBhOfwYAACAYYlaAAAAhiVqAQAAGJaoBQAAYFiiFgAAgGGJWgAAAIYlagEAABiWqAUAAGBY/x9yb0V3A30ZgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_hets_AA_AT = []\n",
    "for sim in range(num_sims):\n",
    "    sim_hets = 0\n",
    "    for ofs in range(20):\n",
    "        sim_hets += 1 if random.choice([0, 1]) == 1 else 0\n",
    "    num_hets_AA_AT.append(sim_hets)\n",
    "    \n",
    "fig, ax = plt.subplots(1,1, figsize=(16,9))\n",
    "ax.hist(num_hets_AA_AT, bins=range(20))\n",
    "print(len([num_hets for num_hets in num_hets_AA_AT if num_hets==20]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd9ae5b79e8>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAIMCAYAAAAq64s8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+wnXV9L/r3pwQltpqARA9NyEl0Qqf+OKKk6h1qx6tHRGVArHrCnBb0eEu1cO/RaWcae2auObZ2co61TO306kCl4h0b/EGRzIhSDnh1Dv4ogaKA6EmwAXZggBKJMEYt+L1/7Cd2GfbO3mTvZO3vzus1s2at9Xm+z9qfNc88O/ud57u+q1prAQAAgB79wrgbAAAAgIMl1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANCtJeNu4GAdf/zxbc2aNeNuAwAAgEPgpptu+ufW2oqZxnUbatesWZNt27aNuw0AAAAOgaq6azbjTD8GAACgW0ItAAAA3RJqAQAA6Fa3n6kFAAA4Ev3Lv/xLJiYm8qMf/WjcrcyLY445JqtWrcrRRx99UPsLtQAAAB2ZmJjI05/+9KxZsyZVNe525qS1loceeigTExNZu3btQb2G6ccAAAAd+dGPfpRnPvOZ3QfaJKmqPPOZz5zTVWehFgAAoDOLIdDuM9f3ItQCAADwpF155ZWpqnznO9/5ufpFF12UY445Jnv27DksffhMLQAAQM8uemGy5+75e71lq5P33DrjsC1btuTXf/3Xc/nll2fTpk0/V/+1X/u1XHnllXnb2942f31NQ6gFAADo2Z67k03zeFV007IZhzz66KO54YYb8qUvfSlnnnnmz0LtnXfemUcffTQf/OAH86d/+qeHJdSafgwAAMCT8rnPfS6nn356TjrppBx33HG5+eabk0xepT3nnHPyile8It/97nfzwAMPHPJehFoAAACelC1btmTDhg1Jkg0bNmTLli1JkssvvzwbNmzIL/zCL+RNb3pTPvOZzxzyXkw/BgAAYNYeeuihXH/99bnttttSVXn88cdTVfmt3/qtbN++Pa95zWuSJD/5yU/ynOc8JxdccMEh7ceVWgAAAGbts5/9bM4999zcdddd2blzZ+65556sXbs27373u7Np06bs3LkzO3fuzL333ptdu3blrrvuOqT9CLUAAADM2pYtW3L22Wf/XO03f/M3s3PnzifUzz777Fx++eWHtJ9qrR3SH3CorF+/vm3btm3cbQAAABxWd9xxR371V3/1Xwtj+kqf+fSE95Skqm5qra2faV+fqQUAAOjZYQ6gC43pxwAAAHRrxlBbVSdW1Zeq6o6qur2q/vNQP66qrq2q7cP9sUO9qurDVbWjqr5VVS8Zea3zhvHbq+q8kfopVXXrsM+Hq6oOxZsFAABgcZnNldrHkvx+a+1Xk7w8yQVV9bwkG5Nc11pbl+S64XmSvC7JuuF2fpKPJJMhOMn7krwsyUuTvG9fEB7GnD+y3+lzf2sAAAAsdjOG2tbafa21m4fHjyS5I8nKJGcluWwYdlmSNw6Pz0ryiTbp60mWV9UJSV6b5NrW2u7W2veTXJvk9GHbM1prX2uTq1Z9YuS1AAAAYFpPaqGoqlqT5MVJvpHk2a21+5LJ4FtVzxqGrUxyz8huE0PtQPWJKeoAP2++V/Y7lMawaiAAwJFo1qG2qn4pyRVJ3t1a+8EBPvY61YZ2EPWpejg/k9OUs3r16plaBhabPXcnm/aMu4vZ2bRs3B0AABwyO3fuzBlnnJHbbrttVuM/97nP5aSTTsrznve8ee9lVqG2qo7OZKD9ZGvt74by/VV1wnCV9oQkDwz1iSQnjuy+Ksm9Q/2V+9X/v6G+aorxT9BauzjJxcnk99TOpncAAIDF7NTN12fXw3vn7fVWLl+aGza+at5eL5kMtWecccZ4Qu2wEvHHktzRWvvzkU1bk5yXZPNwf9VI/cKqujyTi0LtGYLvNUn+dGRxqNOSvLe1truqHqmql2dyWvO5Sf5yHt4bAADAorfr4b3ZufkN8/Z6azZ+flbjHn/88fzO7/xOvvrVr2blypW56qqrcu+99+aCCy7Igw8+mKc97Wm55JJLsnv37mzdujVf/vKX8yd/8ie54oor8tznPnfe+p3NldpTk/x2klur6pah9keZDLOfrqp3JLk7yVuGbVcneX2SHUl+mOTtSTKE1z9OcuMw7v2ttd3D43cl+XiSpUm+MNwAAABYoLZv354tW7bkkksuyVvf+tZcccUV+Zu/+Zt89KMfzbp16/KNb3wjv/d7v5frr78+Z555Zs4444y8+c1vnvc+Zgy1rbX/mak/95okr55ifEtywTSvdWmSS6eob0vygpl6AQAAYGFYu3ZtTj755CTJKaeckp07d+arX/1q3vKWt/xszI9//OND3seTWv0YAAAAkuSpT33qzx4fddRRuf/++7N8+fLccsstB9hr/s34PbUAAAAwk2c84xlZu3ZtPvOZzyRJWmv55je/mSR5+tOfnkceeeSQ/FyhFgAAgHnxyU9+Mh/72Mfyohe9KM9//vNz1VWT6wlv2LAhH/zgB/PiF784d95557z+TNOPAQAAOrZy+dJZr1g829ebyZo1a37uO2r/4A/+4GePv/jFLz5h/Kmnnppvf/vb89PgfoRaAACAjs33d8r2xvRjAAAAuiXUAgAA0C2hFgAAoDOttXG3MG/m+l6EWgAAgI4cc8wxeeihhxZFsG2t5aGHHsoxxxxz0K9hoSgAAICOrFq1KhMTE3nwwQfH3cq8OOaYY7Jq1aqD3l+oBQAA6MjRRx+dtWvXjruNBcP0YwAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANCtJeNuABizi16Y7Ll73F3MzrLV4+4AAIAFRqiFI92eu5NNe8bdBQAAHBTTjwEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3Zgy1VXVpVT1QVbeN1D5VVbcMt51VdctQX1NVe0e2fXRkn1Oq6taq2lFVH66qGurHVdW1VbV9uD/2ULxRAAAAFp/ZXKn9eJLTRwuttf/QWju5tXZykiuS/N3I5jv3bWutvXOk/pEk5ydZN9z2vebGJNe11tYluW54DgAAADOaMdS21r6SZPdU24arrW9NsuVAr1FVJyR5Rmvta621luQTSd44bD4ryWXD48tG6gAAAHBAc/1M7SuS3N9a2z5SW1tV/1hVX66qVwy1lUkmRsZMDLUkeXZr7b4kGe6fNceeAAAAOEIsmeP+5+Tnr9Lel2R1a+2hqjolyeeq6vlJaop925P9YVV1fianMGf16tUH0S4AAACLyUGH2qpakuRNSU7ZV2ut/TjJj4fHN1XVnUlOyuSV2VUju69Kcu/w+P6qOqG1dt8wTfmB6X5ma+3iJBcnyfr16590KAYAft6pm6/Prof3jruNGa1cvjQ3bHzVuNsAYAGay5Xaf5/kO621n00rrqoVSXa31h6vqudkckGo77XWdlfVI1X18iTfSHJukr8cdtua5Lwkm4f7q+bQEwDwJOx6eG92bn7DuNuY0ZqNnx93CwAsULP5Sp8tSb6W5FeqaqKq3jFs2pAnLhD1G0m+VVXfTPLZJO9sre1bZOpdSf46yY4kdyb5wlDfnOQ1VbU9yWuG5wAAADCjGa/UttbOmab+tilqV2TyK36mGr8tyQumqD+U5NUz9QEAAAD7m+vqxwAAADA2Qi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdWjLuBgAWpWWrk03Lxt3FzJatTt5z67i7AAA4aEItwKHQS1DsIXgDAByA6ccAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALq1ZNwNAMBic+rm67Pr4b3jbmNWVi5fOu4WAGBOhFoAmGe7Ht6bnZvfMO42AOCIYPoxAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3Voy7gYAAGaycvnSrNn4+XG3MSsrly/NDRtfNe42AI4YQi0AsOD1FBJ7Cd8Ai4XpxwAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbM4baqrq0qh6oqttGapuqaldV3TLcXj+y7b1VtaOqvltVrx2pnz7UdlTVxpH62qr6RlVtr6pPVdVT5vMNAgAAsHjN5krtx5OcPkX9otbaycPt6iSpqucl2ZDk+cM+/09VHVVVRyX5qySvS/K8JOcMY5Pkvw2vtS7J95O8Yy5vCAAAgCPHjKG2tfaVJLtn+XpnJbm8tfbj1to/JdmR5KXDbUdr7XuttZ8kuTzJWVVVSV6V5LPD/pcleeOTfA8AAAAcoebymdoLq+pbw/TkY4fayiT3jIyZGGrT1Z+Z5OHW2mP71QEAAGBGBxtqP5LkuUlOTnJfkg8N9ZpibDuI+pSq6vyq2lZV2x588MEn1zEAAACLzkGF2tba/a21x1trP01ySSanFyeTV1pPHBm6Ksm9B6j/c5LlVbVkv/p0P/fi1tr61tr6FStWHEzrAAAALCIHFWqr6oSRp2cn2bcy8tYkG6rqqVW1Nsm6JP+Q5MYk64aVjp+SycWktrbWWpIvJXnzsP95Sa46mJ4AAAA48iyZaUBVbUnyyiTHV9VEkvcleWVVnZzJqcI7k/xukrTWbq+qTyf5dpLHklzQWnt8eJ0Lk1yT5Kgkl7bWbh9+xB8mubyq/iTJPyb52Ly9OwAAABa1GUNta+2cKcrTBs/W2geSfGCK+tVJrp6i/r386/RlAAAAmLW5rH4MAAAAYyXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt5aMuwFYtC56YbLn7nF3MbNlq8fdAQAAHDShFg6VPXcnm/aMuwsAAFjUTD8GAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuzRhqq+rSqnqgqm4bqX2wqr5TVd+qqiuravlQX1NVe6vqluH20ZF9TqmqW6tqR1V9uKpqqB9XVddW1fbh/thD8UYBAABYfGZzpfbjSU7fr3Ztkhe01v5dkv+V5L0j2+5srZ083N45Uv9IkvOTrBtu+15zY5LrWmvrklw3PAcAAIAZzRhqW2tfSbJ7v9rft9YeG55+PcmqA71GVZ2Q5Bmtta+11lqSTyR547D5rCSXDY8vG6kDAADAAc3HZ2r/U5IvjDxfW1X/WFVfrqpXDLWVSSZGxkwMtSR5dmvtviQZ7p81Dz0BAABwBFgyl52r6r8keSzJJ4fSfUlWt9YeqqpTknyuqp6fpKbYvR3Ezzs/k1OYs3r16oNrGgAAgEXjoK/UVtV5Sc5I8h+HKcVprf24tfbQ8PimJHcmOSmTV2ZHpyivSnLv8Pj+YXryvmnKD0z3M1trF7fW1rfW1q9YseJgWwcAAGCROKhQW1WnJ/nDJGe21n44Ul9RVUcNj5+TyQWhvjdMK36kql4+rHp8bpKrht22JjlveHzeSB0AAAAOaMbpx1W1JckrkxxfVRNJ3pfJ1Y6fmuTa4Zt5vj6sdPwbSd5fVY8leTzJO1tr+xaZelcmV1JemsnP4O77HO7mJJ+uqnckuTvJW+blnQEAALDozRhqW2vnTFH+2DRjr0hyxTTbtiV5wRT1h5K8eqY+AAAAYH/zsfoxAAAAjIVQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuzbj6MQAsBKduvj67Ht477jZmZeXypeNuAQCOGEItAF3Y9fDe7Nz8hnG3ATNauXxp1mz8/LjbmJWVy5fmho2vGncbAHMi1AIAzKOeQmIv4RvgQHymFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdWjLuBgAYo2Wrk03Lxt3FLP3tuBsAABYgoRbgSPaeW8fdwext/Py4OwAAFiDTjwEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQrVmF2qq6tKoeqKrbRmrHVdW1VbV9uD92qFdVfbiqdlTVt6rqJSP7nDeM315V543UT6mqW4d9PlxVNZ9vEgAAgMVptldqP57k9P1qG5Nc11pbl+S64XmSvC7JuuF2fpKPJJMhOMn7krwsyUuTvG9fEB7GnD+y3/4/CwAAAJ5gVqG2tfaVJLv3K5+V5LLh8WVJ3jhS/0Sb9PUky6vqhCSvTXJta213a+37Sa5Ncvqw7Rmtta+11lqST4y8FgAAAExrLp+pfXZr7b4kGe6fNdRXJrlnZNzEUDtQfWKK+hNU1flVta2qtj344INzaB0AAIDF4FAsFDXV52HbQdSfWGzt4tba+tba+hUrVsyhRQAAABaDuYTa+4epwxnuHxjqE0lOHBm3Ksm9M9RXTVEHAACAA5pLqN2aZN8KxucluWqkfu6wCvLLk+wZpidfk+S0qjp2WCDqtCTXDNseqaqXD6senzvyWgAAADCtJbMZVFVbkrwyyfFVNZHJVYw3J/l0Vb0jyd1J3jIMvzrJ65PsSPLDJG9Pktba7qr64yQ3DuPe31rbt/jUuzK5wvLSJF8YbgAAAHBAswq1rbVzptn06inGtiQXTPM6lya5dIr6tiQvmE0vAAAAsM+hWCgKAAAADguhFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0C2hFgAAgG4JtQAAAHRLqAUAAKBbQi0AAADdEmoBAADollALAABAt4RaAAAAuiXUAgAA0K2DDrVV9StVdcvI7QdV9e6q2lRVu0bqrx/Z571VtaOqvltVrx2pnz7UdlTVxrm+KQAAAI4MSw52x9bad5OcnCRVdVSSXUmuTPL2JBe11v5sdHxVPS/JhiTPT/LLSf5HVZ00bP6rJK9JMpHkxqra2lr79sH2BgAAwJHhoEPtfl6d5M7W2l1VNd2Ys5Jc3lr7cZJ/qqodSV46bNvRWvteklTV5cNYoRYA4BBauXxp1mz8/LjbmNHK5Utzw8ZXjbsNYIGar1C7IcmWkecXVtW5SbYl+f3W2veTrEzy9ZExE0MtSe7Zr/6yeeoLAIBp9BIUewjewPjMeaGoqnpKkjOTfGYofSTJczM5Nfm+JB/aN3SK3dsB6lP9rPOraltVbXvwwQfn1DcAAAD9m4/Vj1+X5ObW2v1J0lq7v7X2eGvtp0kuyb9OMZ5IcuLIfquS3HuA+hO01i5ura1vra1fsWLFPLQOAABAz+Yj1J6TkanHVXXCyLazk9w2PN6aZENVPbWq1iZZl+QfktyYZF1VrR2u+m4YxgIAAMABzekztVX1tEyuWvy7I+X/XlUnZ3IK8c5921prt1fVpzO5ANRjSS5orT0+vM6FSa5JclSSS1trt8+lLwBm59TN12fXw3vH3casrIyPnQAATzSnUNta+2GSZ+5X++0DjP9Akg9MUb86ydVz6QWAJ2/Xw3uzc/Mbxt3G7GxaluRt4+4CAFhg5mP6MQAAAIyFUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbi0ZdwPwpFz0wmTP3ePuYnaWrR53BwAAsOgJtfRlz93Jpj3j7gIAAFggTD8GAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6NaScTcAALOybHWyadm4u5idZauT99w67i4A4Igg1ALQh55CYi/hGwAWAdOPAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuzTnUVtXOqrq1qm6pqm1D7biquraqtg/3xw71qqoPV9WOqvpWVb1k5HXOG8Zvr6rz5toXAAAAi998Xan931trJ7fW1g/PNya5rrW2Lsl1w/MkeV2SdcPt/CQfSSZDcJL3JXlZkpcmed++IAwAAADTOVTTj89Kctnw+LIkbxypf6JN+nqS5VV1QpLXJrm2tba7tfb9JNcmOf0Q9QYAAMAiMR+htiX5+6q6qarOH2rPbq3dlyTD/bOG+sok94zsOzHUpqsDAADAtJbMw2uc2lq7t6qeleTaqvrOAcbWFLV2gPrP7zwZms9PktWrVx9MrwAAACwic75S21q7d7h/IMmVmfxM7P3DtOIM9w8MwyeSnDiy+6ok9x6gvv/Puri1tr61tn7FihVzbR0AAIDOzSnUVtUvVtXT9z1OclqS25JsTbJvBePzklw1PN6a5NxhFeSXJ9kzTE++JslpVXXssEDUaUMNAAAApjXX6cfPTnJlVe17rb9trX2xqm5M8umqekeSu5O8ZRh/dZLXJ9mR5IdJ3p4krbXdVfXHSW4cxr2/tbZ7jr0BAACwyM0p1LbWvpfkRVPUH0ry6inqLckF07zWpUkunUs/AAAAHFkO1Vf6AAAAwCEn1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeWjLsBAAA4kJXLl2bNxs+Pu41ZWbl8aW7Y+KpxtwFHFKEWAIAFraeQ2Ev4hsVEqAU4BE7dfH12Pbx33G3MaOXypeNuAQBgToRagENg18N7s3PzG8bdBgDAomehKAAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0sOdseqOjHJJ5L8myQ/TXJxa+0vqmpTkt9J8uAw9I9aa1cP+7w3yTuSPJ7k/2qtXTPUT0/yF0mOSvLXrbXNB9sXAIzdstXJpmXj7mJ2lq1O3nPruLsAgIN20KE2yWNJfr+1dnNVPT3JTVV17bDtotban40OrqrnJdmQ5PlJfjnJ/6iqk4bNf5XkNUkmktxYVVtba9+eQ28AMD49hcRewjcATOOgQ21r7b4k9w2PH6mqO5KsPMAuZyW5vLX24yT/VFU7krx02Lajtfa9JKmqy4exQi0AAAAHNC+fqa2qNUlenOQbQ+nCqvpWVV1aVccOtZVJ7hnZbWKoTVcHAACAA5pzqK2qX0pyRZJ3t9Z+kOQjSZ6b5ORMXsn90L6hU+zeDlCf6medX1Xbqmrbgw8+ONUQAAAAjiBzCrVVdXQmA+0nW2t/lySttftba4+31n6a5JL86xTjiSQnjuy+Ksm9B6g/QWvt4tba+tba+hUrVsyldQAAABaBgw61VVVJPpbkjtban4/UTxgZdnaS24bHW5NsqKqnVtXaJOuS/EOSG5Osq6q1VfWUTC4mtfVg+wIAAODIMZfVj09N8ttJbq2qW4baHyU5p6pOzuQU4p1JfjdJWmu3V9WnM7kA1GNJLmitPZ4kVXVhkmsy+ZU+l7bWbp9DXwAAABwh5rL68f/M1J+HvfoA+3wgyQemqF99oP0AAABgKvOy+jEAAACMg1ALAABAt4RaAAAAuiXUAgAA0K25rH7MYnHRC5M9d4+7i9lZtnrcHQAAAAuIUMtkoN20Z9xdAAAAPGmmHwMAANAtoRYAAIBNC7dVAAAJP0lEQVRuCbUAAAB0S6gFAACgW0ItAAAA3RJqAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuLRl3AwCzderm67Pr4b3jbmNWVi5fOu4WAACOCEIt0I1dD+/Nzs1vGHcbAAAsIKYfAwAA0C1XagEAYJ6sXL40azZ+ftxtzMrK5Utzw8ZXjbsNmDOhFgAA5klPIbGX8A0zMf0YAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbi0ZdwMAwBgtW51sWjbuLma2bHXynlvH3QUAC5BQCwBHsl6CYg/BG4CxMP0YAACAbgm1AAAAdEuoBQAAoFtCLQAAAN0SagEAAOiWUAsAAEC3hFoAAAC6JdQCAADQLaEWAACAbgm1AAAAdGvJuBsAxuvUzddn18N7x93GrKxcvnTcLQDAorFy+dKs2fj5cbcxKyuXL80NG1817jZYoIRaOMLtenhvdm5+w7jbAAAOs55CYi/hm/Ew/RgAAIBuuVJ7qFz0wmTP3ePuYnaWrR53BwBwYMtWJ5uWjbuL2Vm2OnnPrePuAuCIIdQeKnvuTjbtGXcXALA49BQSewnfAIuE6ccAAAB0S6gFAACgWwsm1FbV6VX13araUVUbx90PAAAAC9+C+ExtVR2V5K+SvCbJRJIbq2pra+3b4+0MAAAYt16+U9f36Y7Hggi1SV6aZEdr7XtJUlWXJzkriVBLt07dfH12Pbx33G3MaOXypeNuAWBxsVIzzLtegmIPwXsxWiihdmWSe0aeTyR52Zh6YQHrJSgmk2Fx5+Y3jLsNAA63nkJiL+EbOtHLFeVkcV1VXiihtqaotScMqjo/yfnD00er6ruHtKu5OT7/tf553E0wK8cnmfdjdVeSeu98v+oR75AcKw4Jx6ofjlUfDt1x+q9T/RnGHDin+nFEH6tO/lb9t7MZtFBC7USSE0eer0py7/6DWmsXJ7n4cDU1F1W1rbW2ftx9MDPHqh+OVT8cq344Vn1wnPrhWPXDsVo8FsrqxzcmWVdVa6vqKUk2JNk65p4AAABY4BbEldrW2mNVdWGSa5IcleTS1trtY24LAACABW5BhNokaa1dneTqcfcxj7qYJk0Sx6onjlU/HKt+OFZ9cJz64Vj1w7FaJKq1J6zHBAAAAF1YKJ+pBQAAgCdNqJ2jqjq9qr5bVTuqauMU259aVZ8atn+jqtYc/i6pqhOr6ktVdUdV3V5V/3mKMa+sqj1Vdctw+7/H0StJVe2sqluH47Btiu1VVR8ezqtvVdVLxtHnka6qfmXkfLmlqn5QVe/eb4zzakyq6tKqeqCqbhupHVdV11bV9uH+2Gn2PW8Ys72qzjt8XR95pjlOH6yq7wy/366squXT7HvA35XMr2mO1aaq2jXyO+710+x7wL8XmV/THKtPjRynnVV1yzT7Oq86ZPrxHFTVUUn+V5LXZPJriW5Mck5r7dsjY34vyb9rrb2zqjYkObu19h/G0vARrKpOSHJCa+3mqnp6kpuSvHG/Y/XKJH/QWjtjTG0yqKqdSda31qb87rjhj4b/M8nrk7wsyV+01l52+Dpkf8Pvw11JXtZau2uk/so4r8aiqn4jyaNJPtFae8FQ++9JdrfWNg9/WB/bWvvD/fY7Lsm2JOsz+Z3xNyU5pbX2/cP6Bo4Q0xyn05JcPyyk+d+SZP/jNIzbmQP8rmR+TXOsNiV5tLX2ZwfYb8a/F5lfUx2r/bZ/KMme1tr7p9i2M86r7rhSOzcvTbKjtfa91tpPklye5Kz9xpyV5LLh8WeTvLqqfMv5YdZau6+1dvPw+JEkdyRZOd6umIOzMvkPVWutfT3J8uE/LhifVye5czTQMl6tta8k2b1fefTfpMuSvHGKXV+b5NrW2u4hyF6b5PRD1ugRbqrj1Fr7+9baY8PTrydZddgb4wmmOadmYzZ/LzKPDnSshr/D35pky2FtikNKqJ2blUnuGXk+kScGpZ+NGf6B2pPkmYelO6Y0TAF/cZJvTLH5f6uqb1bVF6rq+Ye1MUa1JH9fVTdV1flTbJ/NucfhtSHT/4HgvFo4nt1auy+Z/M++JM+aYozza2H5T0m+MM22mX5XcnhcOEwVv3SaKf3OqYXlFUnub61tn2a786pDQu3cTHXFdf/53LMZw2FSVb+U5Iok726t/WC/zTcn+bettRcl+csknzvc/fEzp7bWXpLkdUkuGKYRjXJeLSBV9ZQkZyb5zBSbnVf9cX4tEFX1X5I8luST0wyZ6Xclh95Hkjw3yclJ7kvyoSnGOKcWlnNy4Ku0zqsOCbVzM5HkxJHnq5LcO92YqlqSZFkObuoKc1RVR2cy0H6ytfZ3+29vrf2gtfbo8PjqJEdX1fGHuU2StNbuHe4fSHJlJqdujZrNucfh87okN7fW7t9/g/Nqwbl/31T94f6BKcY4vxaAYYGuM5L8xzbNAiiz+F3JIdZau7+19nhr7adJLsnUx8A5tUAMf4u/KcmnphvjvOqTUDs3NyZZV1VrhysVG5Js3W/M1iT7Vo58cyYXfvC/c4fZ8PmJjyW5o7X259OM+Tf7Pu9cVS/N5Pnx0OHrkiSpql8cFvNKVf1iktOS3LbfsK1Jzq1JL8/kYg/3HeZW+VfT/q+382rBGf036bwkV00x5pokp1XVscNUytOGGodJVZ2e5A+TnNla++E0Y2bzu5JDbL/1HM7O1MdgNn8vcnj8+yTfaa1NTLXRedWvJeNuoGfDqoQXZvIf+6OSXNpau72q3p9kW2ttayaD1P9bVTsyeYV2w/g6PqKdmuS3k9w6soT7HyVZnSSttY9m8j8d3lVVjyXZm2SD/4AYi2cnuXLIQUuS/G1r7YtV9c7kZ8fq6kyufLwjyQ+TvH1MvR7xquppmVzR83dHaqPHynk1JlW1JckrkxxfVRNJ3pdkc5JPV9U7ktyd5C3D2PVJ3tla+z9aa7ur6o8z+Yd4kry/tWaG0SEyzXF6b5KnJrl2+F349eFbFH45yV+31l6faX5XjuEtHDGmOVavrKqTMzmdeGeG34Wjx2q6vxfH8BaOGFMdq9baxzLF+g/Oq8XBV/oAAADQLdOPAQAA6JZQCwAAQLeEWgAAALol1AIAANAtoRYAAIBuCbUAAAB0S6gFAACgW0ItAAAA3fr/Ad387KWfhrrXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x648 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_AAs_AT_AT = []\n",
    "num_hets_AT_AT = []\n",
    "for sim in range(num_sims):\n",
    "    sim_AAs = 0\n",
    "    sim_hets = 0\n",
    "    for ofs in range(20):\n",
    "        derived_cnt = sum(random.choices([0, 1], k=2))\n",
    "        sim_AAs += 1 if derived_cnt == 0 else 0\n",
    "        sim_hets += 1 if derived_cnt == 1 else 0\n",
    "    num_AAs_AT_AT.append(sim_AAs)\n",
    "    num_hets_AT_AT.append(sim_hets)\n",
    "fig, ax = plt.subplots(1,1, figsize=(16,9))\n",
    "ax.hist([num_hets_AT_AT, num_AAs_AT_AT], histtype='step', fill=False, bins=range(20), label=['het', 'AA'])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Balanced output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip\n",
    "import pickle\n",
    "import random\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mendelian_errors = pickle.load(gzip.open('mendelian_errors.pickle.gz', 'rb'))\n",
    "feature_fit = np.load(gzip.open('feature_fit.npy.gz', 'rb'))\n",
    "ordered_features = np.load(open('ordered_features', 'rb'))\n",
    "num_features = len(ordered_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10905732, 541688)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mendelian_errors), len(list(filter(lambda x: x[0] > 0,mendelian_errors.values())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10905732 10364044 541688 4.9670026734564905\n"
     ]
    }
   ],
   "source": [
    "total_observations = len(mendelian_errors)\n",
    "error_observations = len(list(filter(lambda x: x[0] > 0,mendelian_errors.values())))\n",
    "ok_observations = total_observations - error_observations\n",
    "fraction_errors = error_observations/total_observations\n",
    "print (total_observations, ok_observations, error_observations, 100*fraction_errors)\n",
    "del mendelian_errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(541688, 541287)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob_ok_choice = error_observations / ok_observations\n",
    "\n",
    "def accept_entry(row):\n",
    "    if row[-1] == 1:\n",
    "        return True\n",
    "    return random.random() <= prob_ok_choice\n",
    "\n",
    "accept_entry_v = np.vectorize(accept_entry, signature='(i)->()')\n",
    "\n",
    "accepted_entries = accept_entry_v(feature_fit)\n",
    "balanced_fit = feature_fit[accepted_entries]\n",
    "del feature_fit\n",
    "balanced_fit.shape\n",
    "len([x for x in balanced_fit if x[-1] == 1]), len([x for x in balanced_fit if x[-1] == 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(gzip.open('balanced_fit.npy.gz', 'wb'), balanced_fit, allow_pickle=False, fix_imports=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
